#ifndef __LUFactorization__
#define __LUFactorization__
#include<iostream>
#include<cmath>
#include<vector>
using namespace std;

vector<double> LU(vector<double>& d, vector<double>& u, vector<double>& l, vector<double> b) {
	int n = d.size();
	if (u.size() != n-1 || l.size() != n-1) throw "Invalid Size!";
    vector<double> alpha(n);
    vector<double> beta(n-1);
    alpha[0]=d[0];
    beta[0]=u[0]/alpha[0];
	for (int i = 1; i < n-1; ++ i) {
		alpha[i] = d[i] - l[i-1] * beta[i-1];
        beta[i] = u[i] / alpha[i];
	}
    alpha[n-1]=d[n-1] - l[n-2] * beta[n-2];


    vector<double> x(n);
    vector<double> y(n);
    y[0]=b[0]/alpha[0];
	for (int i = 1; i <n; ++ i){
        y[i] = (b[i] - l[i-1] * y[i-1]) / alpha[i];
    }
    x[n-1]=y[n-1];
    for (int i = n-2; i >=0; -- i){
        x[i] = y[i] - beta[i] * x[i+1];
    }
	// for (int i=0;i<alpha.size();i++)
    //     cout<<alpha.at(i)<<endl;	
    // for (int i=0;i<beta.size();i++)
    //     cout<<beta.at(i)<<endl;	
    // for (int i=0;i<y.size();i++)
    //     cout<<y.at(i)<<endl;	
	return x;
}

#endif
